%这个octave脚本可用于简易演示二维波的叠加、衍射、折射、反射、散射与Doppler效应等。
%CC-BY 4

clc
clear

L=5;
dx = 0.01;
dt = 0.001;

% 界面位于 y = 0.6*L
% 介质性质：
% 0: 无,
% 1: 界面两侧波速不同（折射与反射），下侧速度为c(1), 上侧速度为c(2)
% 2: 界面上有一小孔（衍射）
% 3: 界面上有两小孔（双缝干涉）
% 5: 存在大量“杂质”（散射）
interfaceType = 0;

% 波速
c = [2;1];

[x y] = meshgrid(0:dx:L);
n = size(x, 1);

numSource = 1; %波源个数
source(1).T = 0.1; %周期
source(1).x = 0.5*L; %波源位置
source(1).y = 0.5*L;
source(1).vx = 0.6*c(1); %波源速度，用来演示Doppler效应
source(1).vy = 0;

u0 = zeros(n,n);
u1 = zeros(n,n);
u2 = zeros(n,n);

% 设置波速分布矩阵
C  = c(1)*ones(n,n);
if interfaceType == 0
    ;
elseif interfaceType == 1
    C(y>0.6*L) = c(2);
elseif interfaceType == 2
    C(int32(n * 0.6),:) = 0;
    C(int32(n * 0.6),int32(n/2-5):int32(n/2+5)) = c(1);
elseif interfaceType == 3
    C(int32(n * 0.6),:) = 0;
    C(int32(n * 0.6),int32(n*0.45-5):int32(n*0.45+5)) = c(1);
    C(int32(n * 0.6),int32(n*0.55-5):int32(n*0.55+5)) = c(1);
elseif interfaceType == 5
    for i = 1:10
        _x = L*rand();
        _y = L*rand();
        C((x-_x).^2+(y-_y).^2<0.1.^2)=0;
    end
else
    error('Not Implemented')
end
C = C.^2*(dt/dx).^2;

figure();
TICK=1000;
for tick = 1:TICK
    t(tick) = tick*dt;

    diff_i = u1(1:n-2,2:n-1)-2*u1(2:n-1,2:n-1) + u1(3:n,2:n-1);
    diff_j = u1(2:n-1, 1:n-2)-2*u1(2:n-1, 2:n-1) + u1(2:n-1, 3:n);
    u2(2:n-1,2:n-1) =2*u1(2:n-1,2:n-1) - u0(2:n-1,2:n-1) + C(2:n-1,2:n-1).*(diff_i + diff_j);

    % 波源
    for _n = 1:numSource
        _x = uint32((source(_n).x+source(_n).vx*t(tick))/L*n);
        _y = uint32((source(_n).y+source(_n).vy*t(tick))/L*n);

        if(_x < n && _y <n && _x > 0 && _y > 0)
            u2(_y,_x) = sin(2*pi/source(_n).T*t(tick)); %波源
        end
    end

    %固定边界
    u2(1,:) = 0;
    u2(n,:) = 0;
    u2(:,1) = 0;
    u2(:,n) = 0;

	u0 = u1;
	u1 = u2;

    if mod(tick,10) == 0
		clf
		hold on
		axis equal
        axis([0 L 0 L -1 1])
        xlabel('x')
        ylabel('y')
        caxis([-1 1])
        axis off

        view(0,90);
        h = surf(x, y, u2,'EdgeColor','none');
		drawnow
		pause(0.01)
    end
end
